Constraint damping for the Z4c formulation of general relativity 
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One possibility for avoiding constraint violation in numerical relativity simulations adopting free- 
evolution schemes is to modify the continuum evolution equations so that constraint violations 
are damped away. Gundlach et. al. demonstrated that such a scheme damps low-amplitude, high- 
frequency constraint- violating modes exponentially for the Z4 formulation of general relativity. Here 
we analyze the effect of the damping scheme in numerical applications on a conformal decomposition 
of Z4. After reproducing the theoretically predicted damping rates of constraint violations in the 
linear regime, we explore numerical solutions not covered by the theoretical analysis. In particular 
we examine the effect of the damping scheme on low- frequency and on high-amplitude perturbations 
of flat spacetime as well and on the long-term dynamics of puncture and compact star initial data 
in the context of spherical symmetry. We find that the damping scheme is effective provided that 
the constraint violation is resolved on the numerical grid. On grid noise the combination of artificial 
dissipation and damping helps to suppress constraint violations. We find that care must be taken 
in choosing the damping parameter in simulations of puncture black holes. Otherwise the damping 
scheme can cause undesirable growth of the constraints, and even qualitatively incorrect evolutions. 
In the numerical evolution of a compact static star we find that the choice of the damping parameter 
is even more delicate, but may lead to a small decrease of constraint violation. For a large range of 
values it results in unphysical behavior. 



I. INTRODUCTION 



The most common way to construct numerical solu- 
tions to the field equations of general relativity is to take 
a free-evolution approach. The Hamiltonian and mo- 
mentum constraints of the theory are explicitly solved 
only for initial data. Then the remaining field equa- 
tions are rewritten in a suitable hyperbolic form, and 
the initial data can be evolved using the desired numer- 
ical method with this hyperbolic formulation. In the ab- 
sence of boundaries the contracted Bianchi identities can 
be used to show that if the constraints are satisfied on 
one spacelike slice of a foliation, then they will be sat- 
isfied everywhere. However, numerical solutions violate 
the constraints. This violation can be considered under 
control, if when one applies more resolution to the prob- 
lem, the constraint violation converges away at an appro- 
priate rate. Nonetheless, even if the constraint violation 
converges away, at finite resolution constraint violation 
is undesirable. A number of strategies to minimize the 
violation have been considered. One is to choose the 
formulation such that every constraint propagates. In 
combination with suitable boundary conditions, the con- 
straint violation on the numerical grid should then be 
propagated away. On the other hand, if the constraints 
do not propagate then any violation may sit on the grid 
and grow. Another strategy is to use a constraint damp- 
ing scheme, namely to modify the evolution equations 
so that the constraint-satisfying hypersurface becomes 
an attractor in phase space; such an evolution scheme 
is sometimes called a A-systcm [ij. The constraints are 
then referred to as the A-variablcs. 



The Z4 formulation [2H7[ has both propagating con- 
straints and admits a constraint damping scheme [8| . The 
Z4 formulation has a close relationship with the gener- 



alized harmonic formulation, and the damping scheme 
is essentially the same for both systems. The damping 
scheme was a crucial ingredient in the first successful evo- 
lution of orbiting binary black holes through merger [9| . 
Analytic calculations demonstrating that constraint vio- 
lations will be damped away are performed in the frozen 
coefficient approximation. On the basis of these calcula- 
tions one expects that the damping scheme will be effec- 
tive, in numerical applications, on constraint violations 
that are of low amplitude and high frequency in space- 
times that are close to stationarity. Since the constraint 
damping scheme is a modification of the continuum equa- 
tions, this high-frequency should be resolved on the nu- 
merical mesh. It is not clear what effect the damping 
scheme will have on ill-resolved numerical noise. 

A conformal decomposition of Z4, called Z4c, was pro- 
posed [10, [m with the hope of bringing the advantages 
of propagating constraints and the constraint damping 
scheme to the puncture method [l^, [l3| for the evolu- 
tion of binary black holes. Here we continue that in- 
vestigation, considering more carefully the effect of the 
constraint damping scheme on numerical evolutions. We 
address the following questions: (i). Under what condi- 
tions can the theoretically predicted damping rates be 
recovered in the numerical approximation? (ii). How ef- 
fective is the damping scheme in astrophysically relevant 
spacetimes? (iii). In practical applications what are rea- 
sonable values for the constraint damping coefficients? 

A variation of the conformal decomposition has been 
recently presented in |14l |. There it was found that the 
constraint damping terms are essential for stable long- 
term 3D evolutions of binary black holes and the gauge 
wave test. Since the Z4c conformal decomposition differs 
from that in |14l | by nonprincipal terms and implementa- 
tion details (e.g. constraints projection and summation- 



by-part operators), our study does not necessarily apply 
in that case. More work is required to carefully evaluate 
the role of the constraint damping scheme in that case. 

Because of the obvious computational overhead of 
working in three dimensions, we once again present nu- 
merical results in spherical symmetry. Note that since 
the constraint damping scheme is a modification to the 
continuum Z4c formulation, our results arc expected to 
reflect the behavior of the full system. Working in spher- 
ical symmetry furthermore affords us the possibility of 
performing a thorough study in the parameter space of 
constraint damping coefScients. 

In Sec. |IT] we summarize the equations of motion for 
the Z4c formulation and describe the constraint damp- 
ing scheme we employ. We also present the expected 
theoretical rates of damping in the high-frequency, frozen 
coefficient approximation. In Sec. IIIII wc present our nu- 
merical study. Finally wc conclude in Sec. IIVI 

Notation. Geometric units are employed. Standard 
notation for the 3+1 general relativity is used, e.g. 
di, partial derivative with respect to coordinates x', 
i = 1,2,3, Di, 3-covariant derivative, £^, Lie derivative 
along the vector /?*, 7, determinant of the 3- metric 7^, 
Kij, a, P'^ , extrinsic curvature, lapse function and shift 
vector. In the perturbed flat-space simulations there is 
no natural scale of time; there we give the time in abi- 
trary units. 



II. THE Z4C CONSTRAINT DAMPING 
SCHEME 



Their time dependence can be computed as (neglecting 
matter terms). 



dtH = - 2aD'M, - AM.D'a + 2aKH 

-AKi{l + K2)aKQ + Ci3H, (8) 

dtMi = - -aD,H + aKM, - {Dia)H 

+ D^ (2a^(,Z,)) - A (2a7''^(fcZ,)^ 

+ 2Ki(l + K2)A(ae)+/:/3M,. (9) 



From Eq. ([4][5]) one sees that 0, Zi behave as A variables 
if the free parameters ki_2 are properly chosen. Conse- 
quently the Einstein constraint are damped for ki > 
and K2 > —1 0. The constraint subsystem is closed. 
If the constraints Q,Zi and H,Mi are satisfied in one 
hypersurface they will remain satisfied at all times. In- 
troducing the following variables and definitions. 



lij =1 '^lij, X = 7 ^ 

k - 7^^' lU, - 26 , iy = 7~ * (A' 



(10) 



^J ^1^JK), (11) 



2f^Z, + f'l'^'ijkj, Td = l^^^jk, 



(12) 



The Z4c formulation 



the conformal evolutions equations, Z4c, read. 



In 3 -1- 1 form the field equations of the Z4c formulation 
of general relativity for the three-metric and extrinsic 
curvature read 



daij = - 2aK.i.j + Cp-fij, 



(1) 



dtK.j = - D,D,a + a[R,j - 2K,''Kkj + KK,j 

+ 2i)(iZj) - Kl(l + K2)jijQ] 

+ 47ra[7„- (S* - Paum) - 2Sij] + CpKi^ , (2) 
where we use the notation 

D,Zj=l--^lkMr'Z% (3) 

The constraints 0, Zi evolve according to 



dtQ ^a[-H + D'Z, - Ki(2 + K2)e] + £^6, 



(4) 



dtZ, ^a[Ah + D,e - mZ,] + 7^^'at[7"^7y] + P' DjZ„ 

(5) 

where the Hamiltonian and momentum constraints 
H, Mi arc given by 



H =R- KijK'^ 



M,. 



K' - 16^Padm = 0, (6) 

^■'[^^J-7^JA1-87^5, = 0. (7) 



dtX^^x[a{K + 2Q)~D,p\ 



(13) 



da,, = - 2aMj + /3'=7.,,fe + 27fc(./?J;) - fl.jP\ , (14) 
dtk = - D'Dia + a[A,jA'^ + -{K + 29)^] 

-f- 47ra[S' + Padm] + aKi{l - K2)e + P'K^, (15) 

+ a[(/^ + 2e)iy-2i^ifc,] 



2 r 



4- /3^A,,- fe + 2Afc(,/3.';.) - -Ay/3^fc 



dtr 



2A'^a 
1 



,+2a[f;.,i^-^--^i^nn(x),, 
f^2k + 9) J - Snf^S,] + l'^P],k 



2aKi(r-fd* 

^-R~-Ai.A'i 
2 2 ^ 

Stt/sadm - k.i(2 + ^2)9] -I- £/39. 



dtS =a[^R ~ \a^A'' + ^(^^ + 20)' 



(16) 



(17) 
(18) 



Here the intrinsic curvature is written as 



1 



(19) 



4^2—.- 4^2--^'^ 



-^D,xD,x - -^irjD'xDix, (20) 



1 



+ f™(2ff(,f,)fc„ + rLffe,,). (21) 

The equations above are constrained by two algebraic 
expressions, ln(det7) = and j^^Aij = 0, which are 
expHcitly imposed during the numerical evolution. A so- 
lution of the evolution system is a solution of the Ein- 
stein system provided that 0, Zi, H and Mi vanish. In 
our numerical applications we close the system with the 
puncture gauge choice [1^, [l^ , 



dta — — ^LCt K 4- £^a, 



(22) 
(23) 



with /iL = 2/a, /.tg = l/a^ and 77 = 2/Af , where M is 
the ADM mass of the spacetime. Unless stated otherwise, 
we employ the constraint-preserving boundary condition 

of [m. 

For a more thorough introduction to the Z4c formula- 
tion we refer the reader to [lOl . |Tl| . The full formulation 
generically forms a strongly hyperbolic system of partial 
differential equations, except in special cases which we 
do not discuss here. 



B. Theoretical damping rates 

Here we briefly review the results of [g , see also [131 for 
a discussion of stability of the undamped nonlinear con- 
straint system. Consider the subsystem (|3][5l |5]0 when 
the initial data is constraint violating. We consider a 
small constraint-violating perturbation on a background 
which satisfies the Z4c equations of motion. Working in 
the frozen coefficient approximation we can choose coor- 
dinates at a point so that the spatial metric is just that 
of flat-space, and the lapse is unity. The only nontriv- 
ial component of the metric is the shift, which can only 
be made to vanish if we allow ourselves the freedom to 
change the spatial slice, which we take as given ^M,- We 
discard nonprincipal terms involving products of the per- 
turbation and the background curvature, extrinsic curva- 
ture, and matter sources, which is inconsistent with our 
first-order perturbation approach. The inclusion of mat- 
ter sources without background curvature terms should 
give the correct damping rates for test matter fields on 
fiat space. Such an analysis is not expected to give the 
rates for compact stars as evolved in Sec. IIIIl nontriv- 
ial background curvature, extrinsic curvature and matter 



source terms will affect the damping rates, but for a con- 
sistent analysis all of these effects should be considered 
together. Already in the variable coefficient approxima- 
tion such an analysis will be challenging. We will use the 
shift only to illustrate that it does not play any role on 
the damping rates. Besides the inclusion of the shift our 
calculations are exactly the same as those of [g . To start 
with, we make a plane wave ansatz 



e = e 



st-^-iujix"^ 



e. 



st-\riuJiX^ y 



(24) 



for the solution of the constraint subsystem, with com- 
plex s and real w^. We restrict to the special case K2 ~ 0, 
and write ki = k. Rewriting the constraint subsystem 
in fully second order in time form, we find that following 
eigenvalue problem 



A + 2sk + 2iuj/3k iuk 

\ + sk — iujj3k 



e 

Zu, 



0, (25) 



( A + sfc - iujkp ) ( z^ ) = 0, (26) 
must be satisfied, where we write 



A = s^ +uj^ -2iujsP-u^P\ 



(27) 



and Zq, stands for the component of Zi in the w' direc- 
tion, while Za are the components transverse to the unit 
wave vector a)' and /3 = P^Cji . The symbol of this equa- 
tion generically has a complete set of eigenvectors, so the 
system can be rotated to diagonal form, 



Ae \ ( sQ + iLo{Zc, ~ !3Q) 
\zj\ Za, 



A; 



Za 



0, 



0. 



(28) 
(29) 



In one dimension, identifying the direction r, one expects 
that the behavior of the combinations of primitive vari- 
ables. 



u<d = dtO + drZr, 

U^ = Zr 



(30) 
(31) 



will be determined, in the linear regime, by the eigenval- 
ues of the symbol provided that the shift is small. The 
eigenvalues are 



Ae =A + 2sk - 2ikujl3, 
s 



k + iujP ± Vfc^-^, 



(32) 



and 



\z, =A + sk — iujkf], 



s = — — + iujf] ± 




(33) 



and finally 



Xzj. =A + sk — iujkp, 




(34) 



TABLE I: Setting used for the numerical simulations pre- 
sented in this paper. The resolution is given in grid points, 
rout is the coordinate location of the computational bound- 
ary, and QfCFL is the Courant-Friedrichs-Levy factor used in 
the time-stepping. 



Initial data 


Resolution 


rout 


QCFL 


Flat 

Puncture 

Star 


4000 
1000 
2000 


100 a.u. 
50 M 
50 M 


0.5 
0.5 
0.4 



In the low-frcqucncy limit uj <S^ k wc find that 



-2k + iuj/3, 



-k + iuj/3 - 



2 

s ~ iujjS — , (35) 

k 



whereas in the high-frequency limit k <ti uj, we have 



sc^ -fc + iw(^=Fl), sc^--+iuj{/3±l), (36) 

so at lower frequencies half of the modes are damped less. 
In the high-frequency limit, the damping scheme causes 
a exponential decay of the constraints with a decay rate 
oi —k and — fc/2 respectively. 



III. NUMERICAL RESULTS 

In this section we present our numerical results. Spher- 
ical symmetry is assumed. We perform a detailed anal- 
ysis of the damping scheme applied to the evolution of 
flat spacetime with different constraint-violating pertur- 
bations in order to reproduce the analytic results and ex- 
plore the regime not accessible to a pen-and-paper analy- 
sis. We consider then nontrivial initial data composed of 
either punctures or a compact star and evolve them us- 
ing different values for the damping parameters. In this 
case the performance of the damping scheme in a strong- 
field regime on constraint violations related, essentially, 
to different truncation errors are investigated. The code 
employed is described in detail in [M^ . 

Numerical setup. In all numerical simluations we use 
fourth-order finite differences for the discretization of the 
spatial derivatives of the metric fields. For the time inte- 
gration we use Runge-Kutta fourth order in the vacuum 
and puncture tests. In the simulation of the compact 
star we employ the Rungc-Kutta third order in combi- 
nation with a high-resolution shock capturing based on 
the local Lax-Friedrichs flux and the convex-essentially- 
non-oscillatory interpolation for the reconstruction of the 
matter fields. Table |I] summarizes the numerical set- 
tings employed for the results presented; convergence 
tests were run on some cases (see also discussions in the 
following paragraphs). 



A. Perturbed flat space-time experiments 

Initial data and parameter space. Perturbations of 
flat space-time are constructed, depending on which of 
the two eigenmodes, uq, Ui^, we want to analyze, by mod- 
ifying the X variable, 

x(0,r) = l + ^cxp('-^')cos('^rV (37) 

or the F'" variable, 



r''(0,r) = Arexp( -^j cos f — ^ 



(38) 



Simulations employing the first kind of initial data are 
analyzed by looking at the eigenmode Me, while those 
employing the second kind by looking at u„. The eigen- 
modes are Fourier-transformed in space for every time 
step, and their decay is studied by means of the power 
spectral density (PSD) at the frequency '^ = f • 

The parameter space, depicted in Fig. [T] (left panel), 
is spanned by the amplitude A and the frequency i^ of 
the initial constraint violation. We vary also ki = k £ 
[0, 1] but keep for simplicity K2 = 0. The parameter 
6 = 10 (fixed) introduces a length scale to the problem, 
which is useful for tuning the frequency v (number of 
cycles in the period b, sec right panel of Fig. [T]). To 
evaluate the strength of the perturbation, the value of 
the perturbation's amplitude A should be compared with 
unity, see Eq. ([57)) . 

High-frequency, low-amplitude corner. In this region 
of the parameter space the analytical results hold. We set 
J/ = 10 (i> = 1) and A = 10"^, Fig. [2] display the results 
obtained. The numerical data show an exponential decay 
of the PSD at the induced frequency i> for both eigen- 
modes and for different values of the damping parameter 
k. The rate is quantified by linear fitting as displayed 
in Fig. [21 Table [U reports the decay rates for different 
values of k and different grid resolution [n is the num- 
ber of grid points), together with the computed fitting 
error. Almost for every case the analytically predicted 
values, i.e. s sa fc (for uq) and s = f (for u^^), lie within 
the error range of the numerically found number. For 
increasing resolution the error gets smaller. Note that 
the large error is caused by the oscillations of the modes. 
The decay rates of the modes agree much better with the 
analytically predicted values than our conservative error 
estimate suggests (see Fig. [3]) 

From high to low frequency. To explore the low- 
frequency regime, we keep the low amplitude A — 10~^ 
and vary the frequency in the range v = [0, 10]. Analytic 
results indicates that there is no damping for constant in 
space modes (zero- frequency modes). The numerical ex- 
periments show that, for decreasing values of the initial 
perturbation frequency in the range [2, 10], the damping 
scheme remains effective with the analytic exponential 
decay rates. The behavior is displayed in Fig. [31 




R'cquciicy i/ 



No damping 



All modes are damped ex- 
ponentially 




20 30 

r [a.u.] 



FIG. 1: (Left) Parameter space for the constraint violation of flat space-time. The violation can be tuned by its amplitude A 
and frequency v. Analytically known are the low-frequency corner, where the violation is not damped and the high-frequency, 
low amplitude-corner where all modes are damped exponentially. (Right) Shape of the constraint-violating initial data. A 
Gaussian curve with amplitude A and full-width half-maximum b is modulated with the frequency u/b. The parameter u tells 
how many oscillations are within b. The figure shows the constraint violation in x at i = with A = 1, b — 10 and u — 5. 
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FIG. 2: Behavior of the eigenmodes uq (left) and u^i (right) for high-frequency i^ = 
violation. The modes are extracted in the Fourier space at i? = 1. For no damping {k 
k > the modes are damped exponentially with different damping rates. 



10, low-amplitude ^ = 10 constraint 
- 0) the modes are stay constant. For a 



The transition from exponential damping to no damp- 
ing happens after the second octave i/ £ [0,2], Fig. \5\ 
As demonstrated by the plot the transition is smooth 
and quite rapid. The experimental fact, observed here, 
that constraint -violating modes of "almost all" nonzero 
frequencies are killed by the damping scheme can be im- 
portant in numerical relativity simulation. 

From low to high amplitude. Increasing the ampli- 
tude of the perturbation, i.e. moving from a perturbative 
regime to a fully nonlinear situation is a delicate proce- 
dure. Our results can be summarized as follows. 

High-amplitude perturbations, up to A ~ 0.1, are 



damped and the damping rates unaffected. The use of 
progressively higher amplitudes first modifies the damp- 
ing rates (constraint-violating modes are less damped) 
and secondly leads to unphysical results and code fail- 
ures. The maximum amplitude which can be reached 
without changing the damping parameter k depends 
on the frequency of the initial perturbation. For low- 
frequency initial perturbations higher amplitudes arc ef- 
fectively damped by the damping scheme than for high 
frequencies. Furthermore, it is not true that increasing 
k generically allows for higher amplitudes. The depen- 
dence on the damping parameter is not monotonic, the 



X k = 0.25 
+ k = 0.50 
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50 




FIG. 3: Fit of the exponential decay of the eigenmodes us (left) and u^: (right) for high-frequency ;/ = 10 and low-amphtude 
A = 10^'' constraint violation. The analytically predicted decay rates (solid lines) lie in every case within the 68% confidence 
interval of the fits which are represented in the figure by the shaded regions. 
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FIG. 4: The rate of the exponential decay of the eigenmodes ue (left) and u^^ (right) for low-amplitude A — 10 constraint 
violations stays constant in the frequency range v £ [2, 10]. 



optimal value for this problem has been experimentally 
fomid to be fc = 0.5. 

Figure [B] shows that, for fc = 0.5 and i/ = 10, the damp- 
ing rates stay the same as in the low-amplitude case for 
an amplitude range A G [0.0001,0.01]. Setting the am- 
plitude to A = 0.1, the damping is no longer exponential 
and. for higher values, the code gives no reasonable re- 



sult. 

Resolution dependency. Convergence of the results 
has been already reported in Table |ll] and briefly dis- 
cussed; we further comment here focusing on represen- 
tative simulations with k = 0.5, v = 10, A ~ 10"'*' 
with varying resolutions. As shown in Fig. [7] the damp- 
ing effect holds longer for higher resolutions, i.e. if the 
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FIG. 5: In the small frequency-range v £ [0, 2] the transition between exponential damping and no damping happens. The 
figure shows for low-amplitude A = 10~* constraint violations the decay of the eigenmodes u@ (left) and u„ (right). 



TABLE II: Fits results for the decay rates for different resolu- 
tions. The parameters of the initial perturbation are ;/ = 10 
and A = 10"*. 



k 


s {n = 2000) 


s {n = 4000) 


s {n = 8000) 


-^analytic 


0.25 


-0.21 ±0.02 


-0.24 ± 0.02 


-0.25 ±0.01 


-0.25 


0.50 


-0.49 ±0.06 


-0.51 ±0.04 


-0.50 ± 0.04 


-0.50 


0.75 


-0.74 ±0.09 


-0.78 ± 0.08 


-0.77 ±0.05 


-0.75 


1.00 


-0.93 ±0.16 


-1.02 ±0.20 


-1.02 ±0.14 


-1.00 


0.25 


-0.12 ±0.01 


-0.12 ±0.01 


-0.12 ±0.01 


-0.125 


0.50 


-0.24 ±0.02 


-0.24 ± 0.02 


-0.24 ±0.01 


-0.250 


0.75 


-0.35 ±0.04 


-0.35 ± 0.04 


-0.36 ± 0.02 


-0.375 


1.00 


-0.47 ±0.07 


-0.48 ± 0.07 


-0.48 ± 0.03 


-0.500 



frequency of the constraint violation is better resolved, 
then the damping scheme works more effectively. This 
suggests that the damping scheme works as long as the 
frequency of the perturbation is well-resolved and it is 
partially expected since it acts at the continuum level. 
The same effect could have been already anticipated from 
Fig.m which refers to a single resolution but different fre- 
quencies which are resolved differently on the grid. 

Very high frequencies, grid modes and dependency on 
artificial dissipation. The effect of the damping scheme 
on frequencies comparable to or of the order of the nu- 
merical grid (grid modes) is finally studied |3l| . For 
these tests we use k = 0.5 and vary the frequency 
1/ e [10,30] (Note the grid spacing is /i = 0.025.). As 
demonstrated in Fig. |S1 if the frequency v of the pertur- 
bation is increased further to a regime where the signal 
is not well-resolved, the damping becomes progressively 
less effective and deviates from the analytic expectation. 
More importantly, the amount of artificial dissipation |19j 
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FIG. 6: Decay of the eigenmode u^i for high-frequency con- 
straint violation with increasing amplitude A. Between A — 
10~* and A = 10"'^ the damping rate does not change. For 
very high-amplitude A = 10~^ the damping of the eigenmode 
is not exponential anymore and the code does not give phys- 
ically reasonable results. 



plays a significant role. The use of the artificial dissipa- 
tion filters out grid modes and generically attenuates high 
frequencies which are aliased to lower ones. Figure[H](left 
panel) shows that the use of different amounts of dissi- 
pation, (T, quantitatively changes the decay rate of the 
eigenmodes. The higher a is, the higher the damping 
rate. However one cannot expect to use arbitrary large 
values of a, so a balance between k and a and the prop- 
erty of the solution have to be studied case by case by 




X(0, r) = 1 + Acxp ( -— ) rand([-l, 1] 



(39) 
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FIG. 7: The damping effect depends on the resolution of 
the initial constraint violation. The figure shows for high- 
frequency V = 10, low-amplitude A — 10^* the damping of 
the eigenmode ue for different resolutions. For high resolu- 
tion the frequency is well-resolved and therefore the damping 
effects lasts longer than for less well-resolved case. 




Time 



FIG. 8: Increasing the frequency of the initial constraint 
violation with low-amplitude A = 10~* even further to very 
high frequencies. At these high frequencies the constraint 
violation get less resolved which weakens the effect of the 
damping. 



performing convergence tests. In our test case we found 
that a ~ 0.05 roughly reproduce the analytic damping 
rates. 

As a final test we evolve initial data containing random 



Note that, in principle, random noise differs from the 
high-frequency perturbation used before because it has a 
flat spectrum. FigurelH] (right panel) shows that the effec- 
tiveness of the damping scheme depends again strongly 
on the value of a in the artificial dissipation operator 
used. However, differently from the previous case, in this 
case we were not able to recover the analytic damping 
rates 



B. Punctures and compact star experiments 

Single puncture. In order to test the damping ef- 
fect on strong-field evolution, we evolve puncture ini- 
tial data [201 . While the initial data are constraint- 
satisfying, a constraint- violating wave leaving the hole 
and propagating outside is observed in numerical simula- 
tions. This feature is generic and not related to the use 
of Z4, but observed also in BSSNOK evolutions, e.g. [2ll| . 
Note however that the constraint violation is converging 
away with resolution, thus not a continuum feature (the 
constraint subsystem in Z4c does not have superlumi- 
nal speeds). Figure [10] (left panel) shows a snapshot of 
the constraint violation leaving the horizon. During the 
evolution the biggest violation is instead found at the 
puncture, where the solution is not smooth. A priori, 
the frequency of the initial constraint- violating wave, as 
well as the later violation at the puncture, cannot be esti- 
mated, whereas their amplitude is expected to be "small" 
(in the sense that it is converging away). From Fig.fTOlit 
is evident that the frequency of the constraint- violating 
wave spans a certain range of frequencies; in terms of 
the length scale given by the mass, M = 1, we mainly 
observed violation at a peak frequency y ~ 0.5. It can 
be considered as "high" and we expect it to be damped 
since it is within the first octave. 

The numerical evolutions performed with different val- 
ues of k show that the use of the damping scheme gener- 
ically introduces a certain dynamics in the constraints, 
whose values in space oscillates in time around a small 
value close to zero. The evolution of the L2 norm of the 
constraint monitor 



C= ^/IP + WM~+W+Z^, 



(40) 



is reported in Fig. [TT] (left panels) for different values of k. 
In these tests artificial dissipation is used with a — 0.007. 
In all the cases the norm at early times is dominated by 
a violation inside the horizon during the gauge adjust- 
ment which leads to the trumpet solution {22i27i|. The 
initial constraint violation wave is also propagated out 
during this phase, and eventually damped depending on 
the value of k. The bottom left panel shows the norm 
outside the horizon at early times and highlights the ef- 
fect of the damping scheme for several values of fc. At 
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FIG. 9: For initial constraint violations with low amplitude and very high frequency, the artificial dissipation starts to have 
an effect. The high-frequency modes are shifted by the artificial dissipation to lower frequencies which are then damped. The 
figure shows the dependence of the decay of the eigenmode ue on the artificial dissipation parameter a using high-frequency 
i/ = 30 (left) and random noise constraint violation initial data (right). 
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FIG. 10: The puncture and the stuffed puncture initial data violate the constraints. This constraint violations leave the black 
hole horizon. The figure shows the Hamiltonian constraint H at time t = 15M and t — 30M. (Left) Puncture initial data, 
(right) Stuffed puncture initial data. 



later times the amplitude of the oscillations in the con- 
straints amplifies around t ^ 1.25 x 10'^ A/ as shown in 
the top left panel. Depending on the value of k, the am- 
plification is observed to saturate and damp {k < 0.2) 
or to keep on growing, contaminating the numerical so- 
lution (fc > 0.2). In the latter case the code eventually 
fails because the boundary conditions implemented [111 
can not sustain such a large violation. 

To assess the relative importance of boundary con- 
ditions and the constraint damping scheme in our nu- 
merical experiments, we performed long evolutions of 
a single puncture with and without both the damping 
scheme, and either constraint-preserving [ll| or Som- 
mcrfeld boundary conditions. The outer boundary was 



placed at 50 A/. The results are presented in Fig. [T^l 
We find that the use of constraint-preserving boundary 
conditions is more important than that of the damping 
scheme in avoiding violations. Although the use of the 
damping scheme with k = 0.02 reduces the violation by 
a factor of 9 when Sommerfeld conditions are used. This 
test is not expected to be representative of more general 
scenarios in which the outer boundary is placed farther 
out with the same resolution, since then, experimentally 
we find that a smaller constraint violation interacting 
with the Sommerfeld condition results in smaller reflec- 
tions. On the other hand, since the constraint-preserving 
conditions arc found to converge numerically, and the 
Sommerfeld conditions do not, it is expected that at 
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FIG. 11: Time evolution of the norm of the constraint monitor. For the figures on the left puncture initial data was used, for 
the figures on the right stuffed puncture initial data. The long time behavior of both cases is the same which can be seen in the 
pictures on the top. The difference between the two initial data can be seen in the two figures on the bottom which show the 
norm of the constraint monitor outside the horizon for early time. The stuffing introduces a big constraint violation compared 
to the normal puncture data. This violation is damped away. The longtime behavior depends on the value of the damping 
parameter while for damping parameter k < 0.3 the damping scheme does not cause problems in the evolution, damping 
parameter k > 0.3 introduce dynamics to the system which leads to increasing constraint violation and a crash of the code. 



some resolution the constraint violation induced by the 
Sonimerfeld conditions will become dominant, even if the 
outer boundary is placed far out. This has been recently 
pointed out for the case of 3d matter simulations and of 
BSSNOK in Q. 

Stuffed puncture. A second series of tests performed 
is the evolution of puncture initial data stuffed in the 
black- hole interior |2l|, [23, |30| • The hole has been stuffed 
inside the horizon r^ = M/2 at Tox = 0.475 using a 
fourth-order polynomial in the eonformal factor 



ip{0,r) =2.97368 



5.83175 — 



2 / r \4 

- 7.75413 (-) . 

(41) 
The polynomial matches the puncture data at r^x 
up to the second derivatives. The initial data are 
clearly constraint-violating and also show the outgoing 
constraint- violation wave Fig. [10] (right panel). The evo- 
lution of the constraint monitor for different k is reported 



in the right panels of Fig. 1111 Only quantitative differ- 
ences with respect to the puncture case are observed. As 
demonstrated in the right bottom panel of Fig. [11] (note 
the difference in the scale with respect the left panel) , the 
damping scheme is again effective in reducing the outgo- 
ing constraint violation during the initial adjustment. At 
late times the situation is completely analogous to the 
puncture evolution and we do not repeat the description, 
see top right panel. 

Stable star. As a final test, we present a study of 
the influence of the damping scheme on the evolution 
of a stable equilibrium model of a compact star of mass 
A/ = 1.4 Mq described by the ideal gas equation of state. 
Initial data are the same as those employed in previous 
works [ifl, [U and constraint-satisfying. Furthermore 
they provide the exact solution for the evolution prob- 
lem since the system is static. At the typical resolutions 
employed (and without damping, fc = 0) stable evolu- 
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FIG. 12: Long-term evolution of a puncture with the Z4c 
formulation. Different boundary conditions were tested with 
and without damping (with k = 0.02). Sommerfeld boundary 
conditions lead to a high constraint violation, which is effec- 
tively suppressed by the damping scheme. Using Sommerfeld, 
the constraint violation is roughly a factor of 9 lower with and 
without damping at late times. Using constraint-preserving 
boundary conditions, without damping, leads to a late-term 
constraint violation which is 20 times lower than that of un- 
damped Sommerfeld. The damping is less important in the 
case of constraint-preserving boundary conditions. Using the 
standard damping value k = 0.02, the improvement is by only 
a factor of 2. 



tions are obtamed for about 10 ms [t ^ 1400 A/). Trun- 
cation errors trigger radial oscillation in the star, which 
are small amplitude, low frequency in space constraint- 
violating modes, v ^ l/(2i?) where R is the star coordi- 
nate radius. 

In previous investigations on stable stars, which con- 
sidered only one value of the damping parameter, we 
found that the constraint damping terms have a negligi- 
ble effect on the dynamics, while constraint propagation 
made a large difference with the equivalent BSSNOK sim- 
ulations [lO[ . The new results discussed here confirm the 
previous ones for the specific k considered, but they also 
show that for certain values of the damping parameter 
the constraint damping terms are not beneficial in the 
long-term evolutions. 

Figure [T3] shows the evolution of the central rest-mass 
density (left) and of the L2 norm of the constraint moni- 
tor (right) for different values of k employed in the damp- 
ing scheme. For k > 0.3 the constraint damping am- 
plifies the radial oscillations and drives the star to col- 
lapse. Large constraint violations indicate the departure 
from the constraint-satisfying solution space (for clarity 
they arc not shown in the right panel). Smaller values of 
k < 0.06 produce instead an effective constraint damping 
at early times (see right panel Fig.lT^. In the long-term 
however, the evolution without damping scheme (fc = 0) 
is always preferable to the evolutions with k > 0.03. In 
the latter cases a long-term growth is observed, simi- 



lar to that seen in the puncture simulations. The value 
k — 0.02 leads to a constraint damped evolution, but its 
effect is almost negligible. More importantly, this choice 
is not robust in other tests performed. In both simula- 
tions employing a different equation of state, specifically 
a polytrope which produces different truncation errors at 
the surface of the star, and in the migration test of [l3| 
we were not able to identify a value of k which leads to an 
efficient constraint damping. We presented here in detail 
only the test with the most varied outcome. 

While the static stable star is a delicate test (every 
small perturbation causes a departure from the Einstein 
solution) , it provides a specific relevant example in which 
the constraint damping fails for a large choice of damping 
parameters. This indicates that the use of the constraint 
damping scheme without specific investigations is poten- 
tially dangerous. 



IV. CONCLUSION 

In order to expand the body of evidence that a con- 
formal decomposition of the Z4 formulation of general 
relativity [2| may be a useful tool for numerical relativity 
we have presented a detailed study of the effect of the 
constraint damping scheme of Gundlach et al. Q- 

We have attempted to answer three questions, which 
we address here specifically: 

(i). Under what conditions can the theoretically pre- 
dicted damping rates be recovered in the numerical ap- 
proximation ? By studying the evolution of parametrized 
constraint violating-perturbations on top of flat space, 
we first found that the predicted damping rates of Q 
are recovered for well-resolved high-frequency constraint 
violations. Varying the frequency of the constraint vi- 
olation, we found that the analytically predicted expo- 
nential decay is maintained over a large, three-octave, 
range. The cut-off in the effectiveness of the scheme oc- 
curs over a small range at low frequencies. On grid noise, 
unsurprisingly, we find that the predicted damping rates 
are not recovered, although the combination of damping 
and artificial dissipation does help to suppress constraint 
violations. The intuitive explanation for this is that ar- 
tificial dissipation aliases the grid noise to lower frequen- 
cies which are well-resolved, and on which the damping 
scheme is effective. Finally we increased the amplitude 
of the constraint violation. At amplitudes above A ~ 0.1 
the damping scheme becomes increasingly less effective, 
after which numerical integration is often not possible, 
either with or without constraint damping. 

(ii). How effective is the damping scheme in astrophys- 
ically relevant spacetimes? For this part of the inves- 
tigation we began by evolving a single puncture black 
hole. We find that the constraint damping scheme sup- 
presses constraint- violating numerical error leaving the 
black hole horizon, but that it generically introduces a 
dynamical behavior to the constraints. The suppression 
of the violation leaving the horizon is furthermore not 
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FIG. 13: (Left) Time evolution of central rest-mass density for a stable star obtained using different values of k. (Right) Time 
evolution of the constraint monitor corresponding to a subset of the evolutions of the right panel. 



dramatic. For reasonable values of the damping param- 
eter a factor of about 2 or 3 is gained in the norm of 
the constraint violation. If the damping parameters are 
chosen too large, the dynamical behavior induced by the 
scheme causes a large constraint violation to hit the outer 
boundary, which in our tests was placed at 50M, even- 
tually causing a code failure, which we think it may be 
possible to avoid by including constraint dampin g te rms 
in the constraint-preserving boundary conditions |ll| ap- 
propriately. Since mixed puncture-black-hole neutron- 
star initial data are not readily available, "black-hole 
stufhng" has been proposed [2l|, |2^, |30|- Therefore to 
investigate the likely effect of the constraint damping 
scheme in mixed binary evolutions, we evolved a sin- 
gle puncture with a constraint- violating interior. Here 
we find qualitatively the same behavior as in the single 
puncture evolutions, the only difference being that the 
size of the constraint-violating numerical error leaving 
the black- hole is larger. Finally on the question of astro- 
physically relevant spacetimes evolutions of a static star 
were performed. Here we find that using the damping 
scheme is generically of minor benefit and can cause an 
unphysical collapse. 

(iii). In practical applications what are reasonable val- 
ues for the constraint damping coefficients? Our flat- 
space tests demonstrate that higher values of the damp- 
ing parameters are preferable, because then faster rates 
of exponential damping are achieved. On the other hand, 
since our evolutions of compact objects suffer from severe 
problems when the damping parameters are chosen too 
large, we suggest that the damping parameter are chosen 
in the range k £ [0,0.1] for puncture evolutions while for 
matter evolution the safest option to use fc = unless 
specific damping tests are performed. 



In summary, at least for the spherical symmetric sys- 
tems studied within this work, the following statements 
about the constraint damping scheme can be made: Con- 
sidering vacuum spacetimes, the damping scheme may 
be, for carefully chosen damping parameters, a useful tool 
for suppressing constraint violations. This is certainly 
true if there are features in the numerical setup which 
cause large constraint violations, for example, Sommer- 
feld boundary conditions or constraint- violating initial 
data. If no such features are present, the damping scheme 
is not essential and can furthermore affect the physics of 
the system if the damping parameters are taken too large. 
In the evolution of a static compact star our numeri- 
cal evidence indicates that the damping scheme some- 
times leads to a slight decrease of constraint violation. 
On the other hand the damping scheme, in combina- 
tion with some numerical setups, causes growth of the 
constraints; in the special cases we have considered the 
damping scheme is of marginal use. 
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